Granular Gases under Extreme Driving 
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We study inelastic gases in two dimensions using event-driven molecular dynamics simulations. 
Our focus is the nature of the stationary state attained by rare injection of large amounts of energy to 
balance the dissipation due to collisions. We find that under such extreme driving, with the injection 
rate much smaller than the collision rate, the velocity distribution has a power-law high energy tail. 
The numerically measured exponent characterizing this tail is in excellent agreement with predictions 
of kinetic theory over a wide range of system parameters. We conclude that driving by rare but 
powerful energy injection leads to a well-mixed gas and constitutes an alternative mechanism for 
agitating granular matter. In this distinct nonequilibrium steady-state, energy cascades from large 
to small scales. Our simulations also show that when the injection rate is comparable with the 
collision rate, the velocity distribution has a stretched exponential tail. 
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Granular materials are ubiquitous in nature, but nev- 
ertheless, fundamental understanding of the properties of 
granular materials presents many challenges |lH5(- Un- 
derlying these challenges are structural inhomogeneities, 
macroscopic particle size, and energy dissipation, all of 
which arc defining features of granular matter. 

Equilibrium gases have Maxwellian velocity distribu- 
tions. Due to the irreversible nature of the dissipativc 
collisions, granular gases are out of equilibrium. Indeed, 
non-Maxwellian velocity distributions are observed in a 
wide range of experiments in driven gran ular matter in- 
cluding in particular shaken grains j6Hl7j| . In such exper- 
iments, energy is injected over a wide range of scales and 
the measured velocity distribution has a stretched expo- 
nential form. To a large extent, a kinetic theory where 
energy injection through the system boundary is modeled 
by a thermostat successfully describes these nonequilib- 
rium steady-states fl8U2ll ]. 

Furthermore, theoretical studies suggest that the 
steady-state is controlled primarily by the ratio between 
the energy injection rate and the collision rate [HI, [Hj]. 
When the injection rate is much larger than the collision 
rate, the velocity distribution is Maxwellian. However, 
when the injection rate is smaller than the collision rate, 
the velocity distribution is non-Maxwellian, and has a 
stretched-exponcntial tail. 

In this study, we focus on the limiting case where the 
injection rate is vanishingly small and energy is injected 
at extremely large velocity scales [24l [25j . Under such ex- 
treme driving, the injected energy cascades down from 
large velocity scales to small scales and thereby coun- 
ters the dissipation by collisions. Kinetic theory shows 
that in the stationary state, the velocity distribution has 



a power-law high energy tail. These theoretical pre- 
dictions were supplemented by Monte-Carlo simulations 
of the homogeneous Boltzmann equation where spatial 
correlations are ignored. However, such nonequilibrium 
steady-states have yet to be observed using more realistic 
molecular dynamics simulations. In this work we carry 
out extensive molecular dynamics simulations to investi- 
gate the behavior of inelastic hard disks under extreme 
driving. 

The goal of this investigation is to establish whether 
extreme driving is a feasible mechanism for driving gran- 
ular matter. Our main result is that under rare but pow- 
erful injection of energy, a granular gas indeed reaches 
a stationary state that is characterized by a power-law 
velocity distribution. Moreover, our simulations quanti- 
tatively confirm the predictions of the kinetic theory as 
the exponent characterizing the tail of the distribution is 
validated over a wide range of parameters. Our results 
show that extreme driving is a feasible mechanism for 
agitating granular matter. 

Kinetic Theory. Our starting point is the observa- 
tion that the purely-collisional homogeneous Boltzmann 
equation supports stationary solutions (2jJ [27} . The evo- 
lution equation for the velocity distribution /(v) of in- 
clastic hard disks takes the form 



a/(y) 

dt 



dndui du 2 |(ui - u 2 ) ■ n| /(ui)/(u 2 ) 



[<5(v - vi) - <5(v - ui)], 
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with n the impact direction. This Boltzmann equation is 
supplemented by the inelastic collision rule which speci- 
fies the post-collision velocities v 12 as a linear combina- 
tion of the pre-collision velocities Ui,2, 

1 + r 

vi,2=ui, 2 — (ui, 2 - u 2 ,i) • nn. (2) 

In an inelastic collision, the normal component of 
the relative velocity reverses sign and is scaled 



2 



down by the restitution coefficient < r < 1, 
(vi — v 2 ) • n = — r(ui — u 2 ) • n. The energy loss equals 
A£ = -!^|( Ul -u 2 ).n| 2 . 

The collision rule © simplifies to a "fragmentation" 
rule u — > (wi,W2) for collisions involving one extremely 
energetic particle with velocity u and a second, implicit, 
particle with speed much less than |u|. The post-collision 
velocities wi = =^-u • nn and w 2 = u • nn fol- 

low by substituting Ui = and u 2 = u, respectively, 
into ([2]). In the limit |v| — > oo, the nonlinear Boltzmann 
equation |T]) becomes linear and its stationary form is 

0= // dndu|u-n|/(u) [£(v-Wi)+£(v-w 2 )-£(v-u)] 



For arbitrary dimension and for arbitrary collision pa- 
rameters, this linear and homogeneous equation admits 
the power-law solution 



(4) 



■ 1+r ^<T-3 



r(f) 



(5) 



We relate these three quantities by a heuristic argu- 
ment and first note that the energy injection rate is sim- 
ply jV 2 . Also, we anticipate that the velocity distribu- 
tion is truncated at the injection scale V. The energy 
dissipation is dominated by the tail of the distribution 
and is controlled by the upper cut-off V. We estimate 
the dissipation rate T as follows [28j | . 



r ~ p f y v 2 \f (—) vdv~ P V 3 {V/v ) 2 



(6) 



where p is the particle density. In the integrand, the first 
term v accounts for the collision rate, and the second 
-'term v 2 accounts for the energy dissipation in an inelastic 
collision. The integration is performed using the velocity 
distribution Q , and since a < 5 the dissipation is indeed 
dominated by the high velocity tail of the distribution. 
Balancing energy injection with dissipation we obtain a 
relationship between the injection rate 7, the injection 
scale V, and the typical velocity Vq, 



In two-dimensions, the subject of our investigation, the 
exponent a obeys the transcendental equation (2~H 



where 2 _Fi is the hypcrgcomctric function. The exponent 
a grows monotonically with the restitution coefficient. 
The limiting values are a = 4.14922 for completely in- 
clastic collisions (r = 0) and a — 5 in the clastic limit 
(r-H). 

The power-law distribution Q is a stationary solu- 
tion of the linear Boltzmann equation ([3]). Yet, Monte 
Carlo methods show that the full nonlinear Boltzmann 
equation (fT]) docs admit a stationary solution with a tail 
given by ([4]). These numerical solutions are computed 
by injecting energy at a rate that is much smaller than 
the collision rate. In an individual injection event a ran- 
domly chosen particle is given a velocity much larger than 
the typical velocity. Such extreme driving maintains a 
steady-state in which energy injection balances energy 
dissipation. 

The physical mechanism underlying these driven 
steady-states is a cascade in which a high energy parti- 
cle collides with a particle of typical energy yielding two 
high energy particles, each with energies less than that of 
the original high energy particle. These two high energy 
particles produce two more high energy particles, again 
by collisions with the much more numerous particles with 
typical energies, and so on. 

Scaling Analysis. Interestingly, there is a family of 
steady-states generated by extreme driving. If f(v) is 
a stationary solution of ([T]), then Vq 2 J{v/vq) with the 
arbitrary typical velocity Vq is also a stationary solution 
because the collision rule ([2]) is invariant under the scale 
transformation v — > v/vq- The energy injection rate 7, 
the velocity injection scale V, and the typical velocity vq 
are related by the energy balance requirement. 



(7) 



Since the collision rate is proportional to pvQ, the di- 
mensionlcss ratio ip of the injection rate to the collision 
rate scales as a power of the velocity ratio V/vq, 



i> ~ (V/vo) 



3-cr 



(8) 



Since a > 3, we expect a wide power- law range, V 3> vq, 
when the injection rate is much smaller than the collision 
rate, ip <C 1. There is no lower cutoff on the injection rate 
ip, below which power-law is not observed; the smaller is 
ip, the broader the power-law range. When -0 is order 
one, the velocity distribution no longer has a power-law 
tail, and of course, when ip 1, the velocity distribution 
should simply mirror the distribution of injected veloci- 
ties. 

Molecular Dynamics Simulations. We used molecular dy- 
namics [29| to simulate inelastic hard disks in a square 
box with elastic walls. In these event driven simulations, 
upon impact, the velocities of the colliding particles are 
updated according to the collision rule ((2|). Subsequent 
to each collision, wc identify the time and location of 
the next collision. The particles undergo purely ballistic 
motion between two successive collisions. 

We implemented the following velocity-dependent 
restitution coefficient I30l 



h-(l-r)(S n /v c )^ 5 n <v c , 
\r 5 n > v c , 



where 5 n = (vi — V2) • n is the normal component of 
the relative velocity. Here, r is the nominal value of the 
restitution coefficient, valid at large velocities, and v c is 
the cutoff velocity, below which collisions become elas- 
tic. With this realistic restitution coefficient [UHIl, we 
avoid inelastic collapse where an infinite number of colli- 
sions can occur in a finite time (33j . Typically, wc set v c 




FIG. 1: The velocity distribution f(v) versus the velocity v. 
Molecular dynamics simulation results (bullets) are compared 
with the power-law tail predicted by kinetic theory (solid 
line) . 
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FIG. 2: The exponent a versus the restitution coefficient r. 
The molecular dynamics results (bullets) are compared with 
the kinetic theory predictions (solid line). 



much smaller than the typical velocity vo, but for small 
restitution coefficients, we must set v c comparable to vq 
to avoid inelastic collapse. 

To maintain a steady-state, we periodically boost a sin- 
gle randomly-selected particle to a large, random veloc- 
ity. These injection events are rare and they are governed 
by a Poisson process with rate 7, that is, with probabil- 
ity 7 dt injection is implemented during the time interval 
[t,t + dt]. The injection speed is selected from a Gaus- 
sian distribution with zero mean and standard deviation 
V. By taking a long-time average, we confirmed that 
the total energy approaches a constant, and hence, that 
the system reaches a statistical steady-state where energy 
injection and energy dissipation balance. Moreover, the 
velocity distributions were produced by sampling particle 
velocities at a very large number (10 8 ) of equally spaced 
time intervals. We stress that the velocities are sampled 
at time intervals that are completely uncorrelated with 
either collision events or injection events. We tested that 
our sampling produces robust velocity distributions, and 
that the velocity distribution, representing an average 
over the entire system, is truly stationary. In particular, 
the system does not enter the homogeneous cooling state 
34j in between the rare injection events. 

We performed numerical simulations using a system of 
N = 10 3 identical particles with diameter 2R = 1 in a 
square box of size L = 400, corresponding to the low area 
fraction = Nn (R/L) 2 = 4.9 x 10~ 3 . Unless noted oth- 
erwise, these parameters are used throughout this study. 
Energy was injected at rate 7 = 5 x 10™ 7 and the injection 
scale was V = 850. First, we considered weakly-inelastic 
particles, r = 0.9, with the cutoff v c = 0.1. With these 
parameters, the injection rate is much smaller than the 
measured collision rate and their ratio is ip = 1.5 x 10 -5 . 
Consequently, there is substantial scale separation be- 
tween the typical velocity vq and the injection velocity 
V. Over this range, the steady-state velocity distribu- 
tion obtained by the molecular dynamics simulations has 
a power-law high energy tail as in Q and the exponent 



a = 4.74 is in very good agreement with the kinetic the- 
ory prediction given by (0, a = 4.74104 (see figure H}. 
We also confirmed that the ratio V/vo « 10 2 is consistent 
with the scaling estimate (J5J. 

Next, we varied the restitution coefficient and re- 
peated the simulations. Over the entire parameter range 
0.1 < r < 0.9, we find stationary velocity distributions 
with a power-law tail. In general, the ratio V/vq is con- 
sistent with the scaling relation ([5]). Moreover, the ex- 
ponent a obtained from the molecular dynamics simu- 
lations is in excellent agreement with the kinetic theory 
predictions (JSJ for all restitution coefficients (figure 0) . 
We thus arrive at our main result that at least for dilute 
gases, extreme driving in the form of rare but powerful 
energy injection generates a steady-state with a broad 
distribution of velocities. The tail of the velocity distri- 
bution is power-law and the characteristic exponent is 
nonunivcrsal as it depends on the restitution coefficient. 

The fact that kinetic theory holds shows that, to good 
approximation, the gas is well-mixed. We comment that 
it is remarkable that extreme driving results in a well- 
mixed gas. On short time scales, energy injection clearly 
generates spatial correlations because a just-energized 
particle transfers much of its energy to nearby particles 
by inelastic collisions. Yet, on larger time scales energetic 
particles break coherent structures which are known to be 
the consequence of inelastic collisions [35j | . While these 
two mechanisms have opposite effects, the simulations in- 
dicate that when a long time average is taken, the latter 
effect dominates. Thus, energy injected at extreme ve- 
locity scales in a tightly localized region of space, ends 
up evenly distributed throughout the system. 

Snapshots of the time evolution of the system following 
energy injection demonstrate how the inelastic cascade 
works (figure [3]). In the initial stages of the cascade, 
the injection affects only a small region in space and, 
moreover, there are strong spatial correlations between 
the velocities of the particles (figure |3J a-c). However, 
after many inelastic collisions, the injected energy ends 
up evenly distributed throughout the system (figure [3] 




FIG. 3: The inelastic energy cascade. Shown are four time- 
ordered snapshots of the gas shortly after an injection event. 
The top three figures show a small window around the injec- 
tion event in the early stages of the cascade: (a) the initial 
energetic particle (red online) (b) two energetic particles af- 
ter one collision, (c) four energetic particles after three col- 
lisions. Figure (d) shows the entire system, with energetic 
particles shown larger (red online), in a late stage of the cas- 
cade. After many collisions, the injected energy is evenly dis- 
tributed throughout the system. The simulation parameters 
are r = 0.8, N = 1000, L = 400, V = 707 and 7 = 2 x 10" 6 . 



d). When a long time average is taken over many en- 
ergy injection events at different locations, the system 
is maintained in a homogeneous, well-mixed state. Spa- 
tial correlations induced by inelasticity and the injection 
mechanism do not affect the predicted power-law velocity 
distributions. 

Therefore, rare, powerful, and spatially localized en- 
ergy injection is a unique mechanism of agitating gran- 
ular gases. This mechanism induces an extended energy 
cascade which distributes the injected energy to the rest 
of the system. This physical mechanism is different than 
energy injection by walls 0, @] or by an effective ther- 
mostat [20j or by multiplicative driving Q where the in- 
jected energy directly affects only a small region of space. 
In this sense, the energy cascade represents a novel mech- 
anism for agitating granular matter. 
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FIG. 4: The velocity distribution at three different densities 
(solid lines). The simulation parameters are: r = 0.8, 7 = 
2 x 10~ 6 , and V = 707. Eq. © with a = 4.57246 is also 
shown as a reference (dashed line). A best fit to a power- law 
yields a = 4.6, 4.6, and 4.7 for cj> = 7.7 x 10~ 5 , 4.9 x 10~ 3 , 
and 7.9 x 10 -2 , respectively. 



Figure [3] also illustrates that the velocity distribution 
is correlated with injection times for our system so that 
the predicted power-law distribution arises only after av- 
eraging over many measurements taken at times that are 
uncorrelated with the injection times. However, for very 
large systems driven with a fixed but small injection rate 
per particle there would be many temporally overlapping 
but spatially well-separated injection events. At any in- 
stant of time in a very large system cascades at all stages 
of development would be present somewhere in the sys- 
tem and the power law tail would be time-independent. 

We performed additional simulations to test whether 
the results arc- robust with respect to change of parame- 
ters. In particular, we varied the area fraction by fixing 
the number of particles and varying the system size. The 
results shown in figure 0] arc for three different area frac- 
tions: <p = 7.7 xlCT 5 , (j) = 4.9xlCT 3 , and0 = 7.9 xlCT 2 . 
The corresponding values of ip are 4.7 x 10~ 4 , 5.4 x 1CP 5 , 
and 3.4 x 10~ 6 , respectively. Note that in all cases ^<1. 
We find the same power-law tail in all three cases, and 
the exponent is in good quantitative agreement with the 
kinetic theory prediction (figure|4|). Thus, the energy cas- 
cade mechanism can be realized even at area fractions as 
high as <p> « 1CP 1 and with ip as small as order 10 -6 . 

We also studied the dependence on the ratio ip be- 
tween the injection rate and the collision rate by varying 
the injection rate 7 and the injection velocity V. In ac- 
cord with ([5]), we find that the range [vq, V] of power-law 
behavior shrinks as ip increases. As long as ip is suffi- 
ciently small, the distribution has a power-law tail (figure 
[5J. When this ratio becomes sufficiently large, the tail 
is no longer algebraic and a sharper decay occurs. For 
ip = 5.8 x 10~ 2 , we find the stretched exponential tail 
f(v x ) ~ exp ( — const, x \v x \^) with ( = 1.52 (figured]). 
This value is consistent with the theoretical value £ = 3/2 
for inelastic gases driven by white noise [2(| Hl| and the 
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FIG. 5: The velocity distribution f(w) versus the normalized 
velocity variable w = v/vo with different ratios of injection 
rate to collision rate: ip = 1.4 X 1(T 3 (bullets), 2.7 x 1(T 3 
(squares), 1.1 x 10 -2 (diamonds), 2.3 x 10 -2 (up-triangles), 
and 4.4 x 10~ 2 (down-triangles). The restitution coefficient 
is r — 0.8. Also shown is the reference theoretical curve (Q 
with a = 4.57246 (dashed line). 
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FIG. 6: The distribution f(v x ) of the horizontal component 
of the velocity v x at a moderate ip — 5.8 x 1CP 2 . These 
simulations are performed with r = 0.8, 7 = 6 x 10~ 4 , and 
V — 1.41. The solid line is a best-fit to the stretched expo- 
nential f(v x ) ~ exp ( — const, x \v x \ ) with £ = 1.52. 



experimental value observed in vigorously shaken beads 
10]. Indeed, in this intermediate injection rate regime, 
the energy cascade becomes localized, and frequent, small 
injections are similar to white noise driving. On the other 
hand, these stretched exponential tails do not relate to 
those observed in [25| as the system cools down after in- 
jection is turned-off. We find stretched exponential tails 
for the range 10 _1 J$ ip j§ 1. When the injection rate 
exceeds the collision rate (ip 3> 1) the entire distribution 
becomes Maxwellian as the velocity distribution simply 
mirrors the distribution of injected velocities. 

Finally, we mention that we even varied the energy 
injection mechanism itself. In particular, to implement 
injection strictly at large energy scales, we used an "en- 
ergy loss counter" to keep track of the total energy dis- 
sipated by collisions since the last injection. When the 
dissipated energy equals a fixed large value, we inject 
this amount of energy into a single randomly chosen par- 
ticle [24], [H|. Using this variant of extreme driving to 
maintain the steady-state, we found similar power-law 
velocity distributions. From these studies, we conclude 
that the parameter ip controls the velocity distribution, 
and that different energy injection mechanisms lead to a 
power-law distribution with the very same exponent a, 
as long as ip <C 1. 

Conclusions. Extensive event-driven simulations show 
that under extreme driving in the form of rare but pow- 
erful energy injections, an inelastic gas reaches a steady- 
state with a broad distribution of energies. Such driven 
steady-states are observed for a wide range of collision 
parameters, densities, and energy injection rates, as long 
as the injection rate is much smaller than the collision 
rate. The velocity distributions have a power-law tail and 
the characteristic exponent is in good agreement with the 
kinetic theory predictions. When the ratio between the 
energy injection rate and the collision rate becomes suffi- 



ciently large, the velocity distribution has a much sharper 
stretched exponential tail. 

We conclude that extreme driving where energy is in- 
jected only at very large scales presents an alternative 
mechanism for agitating granular matter and that such 
driving leads to a fundamentally different steady-state 
compared with traditional driving where energy is in- 
jected over all scales. Realizing this driving in experi- 
ments is a challenge because the agitation must be ap- 
plied only at very large velocities. One possible mech- 
anism is shooting very fast particles into the system. 
While such a system would involve a growing number of 
particles, the injection of energetic particles should lead 
to transfer of energy from large scales to small scales by 
a cascade of collisions. 



*** 



FIG. 7: The time development, left to right, of a cascade in a 
denser system (cj> = 7.9 x 10 -2 and the same parameters as in 
Fig- SI) Particles with speeds greater than a fixed threshold 
are shown larger (red online). 



Interestingly, the power-law velocity distribution ap- 
pears to hold in systems that include dense clusters. We 
have observed that when the density is increased, there is 
a tendency for clustering near the walls with no measur- 
able deviation from the predicted power law (figure [7]) . 
In this case injection leads to explosive breakup of dense 
regions. It is intriguing that the power-law tail is quite 
robust, and extends to situations where the assumptions 
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underlying the kinetic theory approach can no longer be 
justified. 

In our simulations, inelastic collapse does not play a 
role because collisions become elastic at small relative 
velocities. Yet, if the collisions are purely inelastic, there 
should be a competition between the formation of high- 
density regions by inelastic collisions and the destruction 
of such clusters by high energy particles. Elucidating 
this competition is another possibility for further inves- 



tigation. Nonetheless, our simulations suggest that ex- 
treme driving generates well-mixed steady-states despite 
the fact that this driving is very inhomogeneous in both 
space and time. 
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